import math

def pca(data, n_components):
    """
    主成分分析(PCA)实现
    参数:
        data: 输入数据，形状为(n_samples, n_features)
        n_components: 保留的主成分数量
    返回:
        reduced_data: 降维后的数据，形状为(n_samples, n_components)
        explained_variance_ratio: 各主成分解释的方差比例
    """
    # 步骤1: 数据中心化（去均值）
    n_samples = len(data)
    n_features = len(data[0])
    
    # 计算每个特征的均值
    mean = [sum(col) / n_samples for col in zip(*data)]
    
    # 中心化数据
    centered_data = [
        [data[i][j] - mean[j] for j in range(n_features)]
        for i in range(n_samples)
    ]
    
    # 步骤2: 计算协方差矩阵
    cov_matrix = [[0.0 for _ in range(n_features)] for _ in range(n_features)]
    for i in range(n_features):
        for j in range(i, n_features):  # 利用对称性计算上三角，再复制到下三角
            cov = sum(centered_data[k][i] * centered_data[k][j] for k in range(n_samples)) / (n_samples - 1)
            cov_matrix[i][j] = cov
            cov_matrix[j][i] = cov
    
    # 步骤3: 计算协方差矩阵的特征值和特征向量
    eigenvalues, eigenvectors = jacobi_eigenvalue(cov_matrix, eps=1e-10)
    
    # 步骤4: 按特征值降序排序特征向量
    sorted_indices = sorted(range(len(eigenvalues)), key=lambda k: eigenvalues[k], reverse=True)
    sorted_eigenvectors = [eigenvectors[i] for i in sorted_indices]
    sorted_eigenvalues = [eigenvalues[i] for i in sorted_indices]
    
    # 步骤5: 选择前n_components个特征向量构建投影矩阵
    projection_matrix = [
        [sorted_eigenvectors[i][j] for i in range(n_components)]
        for j in range(n_features)
    ]
    
    # 步骤6: 将数据投影到新的特征空间
    reduced_data = [
        [sum(centered_data[i][j] * projection_matrix[j][k] for j in range(n_features)) 
         for k in range(n_components)]
        for i in range(n_samples)
    ]
    
    # 计算解释方差比例
    total_variance = sum(sorted_eigenvalues)
    explained_variance_ratio = [ev / total_variance for ev in sorted_eigenvalues[:n_components]]
    
    return reduced_data, explained_variance_ratio


def jacobi_eigenvalue(matrix, eps=1e-10, max_iter=1000):
    """雅可比迭代法求解对称矩阵的特征值和特征向量"""
    n = len(matrix)
    eigenvectors = [[0.0]*n for _ in range(n)]
    for i in range(n):
        eigenvectors[i][i] = 1.0  # 初始化为单位矩阵
    
    for _ in range(max_iter):
        # 找到最大非对角线元素
        max_val = 0.0
        p, q = 0, 0
        for i in range(n):
            for j in range(i+1, n):
                if abs(matrix[i][j]) > max_val:
                    max_val = abs(matrix[i][j])
                    p, q = i, j
        
        if max_val < eps:
            break  # 收敛条件
        
        # 计算旋转角度
        theta = 0.5 * math.atan2(2 * matrix[p][q], matrix[q][q] - matrix[p][p])
        cos_t = math.cos(theta)
        sin_t = math.sin(theta)
        
        # 旋转矩阵作用于原矩阵
        for i in range(n):
            if i != p and i != q:
                a_ip = matrix[i][p]
                a_iq = matrix[i][q]
                matrix[i][p] = a_ip * cos_t - a_iq * sin_t
                matrix[p][i] = matrix[i][p]
                matrix[i][q] = a_ip * sin_t + a_iq * cos_t
                matrix[q][i] = matrix[i][q]
        
        # 更新p和q行/列
        a_pp = matrix[p][p]
        a_qq = matrix[q][q]
        a_pq = matrix[p][q]
        matrix[p][p] = a_pp * cos_t**2 + a_qq * sin_t**2 - 2 * a_pq * sin_t * cos_t
        matrix[q][q] = a_pp * sin_t**2 + a_qq * cos_t**2 + 2 * a_pq * sin_t * cos_t
        matrix[p][q] = matrix[q][p] = 0.0
        
        # 更新特征向量
        for i in range(n):
            v_ip = eigenvectors[i][p]
            v_iq = eigenvectors[i][q]
            eigenvectors[i][p] = v_ip * cos_t - v_iq * sin_t
            eigenvectors[i][q] = v_ip * sin_t + v_iq * cos_t
    
    # 提取特征值（对角线元素）
    eigenvalues = [matrix[i][i] for i in range(n)]
    return eigenvalues, eigenvectors


# 示例用法
if __name__ == "__main__":
    # 生成示例数据 (4个样本，3个特征)
    data = [
        [1, 2, 3],
        [4, 5, 6],
        [7, 8, 9],
        [10, 11, 12]
    ]
    
    # 降维到2个主成分
    reduced_data, evr = pca(data, n_components=2)
    
    print("降维后的数据:")
    for row in reduced_data:
        print([round(x, 4) for x in row])
    
    print("\n解释方差比例:", [round(x, 4) for x in evr])
    print("累计解释方差比例:", round(sum(evr), 4))